A Continuum Theory for Unstructured Mesh 
Generation in Two Dimensions 



Guy Bunin 

Department of Physics, Technion, 
Haifa 32000, Israel 
buning@tx. technion. ac. il 



Abstract 

A continuum description of unstructured meshes in two dimensions, both for planar 
and curved surface domains, is proposed. The meshes described are those which, 
in the limit of an increasingly finer mesh (smaller cells), and away from irregu- 
lar vertices, have ideally-shaped cells (squares or equilateral triangles), and can 
therefore be completely described by two local properties: local cell size and local 
edge directions. The connection between the two properties is derived by defining 
a Riemannian manifold whose geodesies trace the edges of the mesh. A function <f>, 
proportional to the logarithm of the cell size, is shown to obey the Poisson equation, 
with localized charges corresponding to irregular vertices. The problem of finding a 
suitable manifold for a given domain is thus shown to exactly reduce to an Inverse 
Poisson problem on </>, of finding a distribution of localized charges adhering to the 
conditions derived for boundary alignment. Possible applications to mesh generation 
are discussed. 
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1 Overview 



A mesh is a partition of a domain into smaller parts, typically with simpler 
geometry, called cells. In two dimensions, both on the plane and on curved 
surfaces, cells are usually triangles or quadrilaterals. The shapes of the cells 
may be important; for many applications, cells with shapes similar to an equi- 
lateral triangle or a square are preferred. The problem of mesh generation 
can then be seen as an optimization problem: to find a partition of a domain 
into well-shaped cells, possibly under additional demands, such as cell size 
requirements. 
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Fig. 1. Unstructured vs. structured meshes, (i) Input domain, (ii) Unstructured 
mesh. Irregular vertices are marked, (hi) Structured mesh created by mapping a 
regular grid. 

The mesh generation problem has been the subject of extensive research. Many 
techniques for creating meshes exist, especially in two dimensions [1],[2],[3]. 
Nevertheless, some of the most popular techniques, which create good meshes 
in many cases, are heuristic in nature, and may create less than optimal meshes 
for some inputs. Two of the inherent characteristics of the mesh generation 
problem seem to make it very difficult to solve: 

(Characteristic 1) The constraints on cells' shapes are global. That is, the 
shape of one cell is constrained by the possible shapes of its neighbors, which 
in turn are constrained by the shapes of their neighbors, and so on. Thus, at 
least in principle, the constraints on the mesh layout extend over the whole 
domain (or more precisely, over each connected component of the domain) . 

(Characteristic 2) The problem combines discrete and continuous aspects. 
The number of cells and the mesh connectivity (i.e.: which cells are neigh- 
bors? which faces do they share?) are of discrete nature, whereas the lo- 
cations of the vertices can vary continuously. These aspects are closely in- 
tertwined, preventing the sole use of purely discrete techniques (algebraic, 
graph-theoretic, etc.), or techniques designed for use in problems of contin- 
uous nature. 

The meshes created by mesh generation algorithms can be divided into struc- 
tured meshes, and unstructured meshes. A structured mesh is a mesh whose 
connectivity is that of a regular grid, see Fig. 1 ,(ih). By assuming the con- 
nectivity of the mesh beforehand, the problem of creating such a mesh is 
considerably simplified (see Characteristic 2 above), and reduces to the prob- 
lem of assigning locations to the mesh vertices. One way of doing that is by 
finding a mapping function that maps the domain of the regular grid to the 
domain to be meshed, and using it to map the vertices. Such techniques are 
known as mapping techniques. For small enough cells, the differential prop- 
erties of the mapping function at the cell's location dictate its shape. If, for 
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example, a mapping is angle preserving, then the inner angle of a small enough 
cell will be approximately preserved under the mapping. For a survey, see Ref. 
2. Mapping techniques have also been offered for unstructured meshes, (in the 
context of smoothing a given mesh see [5], [6] and references therein), as long 
as the connectivity of the mesh is given beforehand. 

Just as a structured mesh can be imagined as created by mapping of a region 
of the plane to the domain to be mapped, an unstructured mesh in two di- 
mensions can be imagined as a surface, that is mapped onto the domain to be 
meshed. The simplest example is that of an unstructured mesh with just one 
irregular vertex (a vertex that has more or less than four cells incident upon 
it). The surface to be mapped in this case is a cone. This can be visualized 
using the Volterra construction [7] : Consider a piece of paper with an angular 
section cut out, see Fig. 2,(i). If the two edges of the section are identified, i.e. 
glued together, the paper will assume the shape of a cone, see Fig. 2,(ii). If 
the cone is then mapped to the plane, a regular grid drawn on the cone would 
be mapped to an unstructured mesh such as the one shown in Fig. 2,(iii). The 
mapping shown in Fig. 2,(iii) has the special property of being conformal: a 
small square on the cone is approximately mapped to a square on plane. This 
creates well shaped cells in the resulting mesh: in the limit of an increasingly 
smaller cell, its square shape is preserved under the conformal mapping. A 
similar construction can be imagined for creating an unstructured mesh with 
more than one irregular vertex; each irregular vertex will then correspond to 
one "cone tip" of the surface. 

The approach of the present work to the problem of creating unstructured 
meshes can be expressed as follows: given a domain to be meshed, what surface, 
with a square grid drawn upon it, can be mapped conformally into this domain? 
Thus, it is not just the mapping function that is sought after, but rather the 
surface to be mapped together with the mapping function. Unlike mapping 
techniques, however, the mapping function is required to be conformal. 
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Fig. 2. Cone point, (i) The Volterra construction, (ii) A cone created by the Volterra 
construction, (iii) The cone mapped onto the plane. Alternatively: geodesies on a 
manifold containing a conical singularity. 
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A mathematical framework highly suitable for dealing with such questions is 
Riemannian geometry. Given the domain to be meshed, Riemannian geome- 
try allows one to define a "new geometry" for that domain. This includes a 
redefinition of the distances between points of the domain, and it is used here 
to redefine distances such that a cell edge have unit length. Thus, instead of 
defining the surface that is mapped and the mapping function separately, both 
are treated together, as the mapping induces a new distance definition on the 
domain to be meshed. For example, the cone in Fig. 2,(ii) and the surface in 
Fig. 2 ,(iii) have the "same geometry" (i.e. are isometric) if the distances in 
Fig. 2,(iii) are defined such that cell-edges have unit length. Using the termi- 
nology of Riemannian geometry, the problem can be restated as assigning a 
new metric to the domain being meshed, having the following properties: 

(Property 1) The metric is locally flat everywhere, except at some points, 
called cone points. A "flat" region can be imagined as a bent, but not 
stretched, piece of paper. The cone points are the "tips of the cones" as 
described above. 

(Property 2) A single real function <p is defined in the domain (except at 
cone points), such that at any given point p, the new metric cjij (p) at p 
is proportional to the original metric g%j (p): g%j (p) = e 2 ^gij(p), with 
1 < z, j < 2,. (On the plane with Cartesian coordinates g^- = S^.) We stress 
that the proportionality e 2 ^ can vary throughout the domain. The quantity 
e~^ is the local change in the distance definition and will be associated 
with the local cell size. (In manifold theory terminology, the metric g^ is 
conformally related to gij). 

These two properties are not sufficient for our purposes. First of all, as the 
Volterra construction implies, for a grid of squares to be drawn around a 
cone tip, the total angle as measured around the tip must be a multiple of 
7r/2 radians. Secondly, the mesh must be aligned along the boundary, see 
Fig. 3. To formalize these demands, we define the direction of mesh-edges at 
each point on the surface. Since the edges are assumed to form right angles at 
incidence, this direction is defined up to an addition of ir/2 radiano. The four 
directions at each point will hence be called a cross, and the field of directions 
on the entire domain being meshed will be called a cross-field. The cross-field is 
related to the new metric by requiring that the curves, generated by following 
the directions of the cross-field, along which the edges will be laid, will be 
geodesies of the new metric fy. (Geodesies are the generalizations of straight 
lines for surfaces and manifolds.) The boundary alignment is formalized by 
requiring that the crosses be aligned with the boundary. 

We therefore add the following property to the required properties of the new 



In the case of a triangular mesh, the direction will be defined up to an addition 
of vr/3. 
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metric: 



(Property 3) A cross-field exists. (Exact definition is given in section 3.) 

Property 1 states that the Gaussian curvature K of the manifold with metric 
Qij be identically zero everywhere, except at cone points. Combined with the 
equation g^ = e 2 ^gij of Property 2, the two formulas give a remarkably simple 
result, viz. that must obey the Poisson equation V 2 = K everywhere exept 
at cone points, where K is the Gaussian curvature of the surface. This is a 
well known result in conformal geometry, see section 3. 



The main result of the present work is that for a cross-field to exist the func- 
tion 0, which uniquely defines the manifold, has to obey the Poisson equation 
V 2 = K + p, with K the Gaussian curvature of the surface to be meshed 
(K vanishes for the plane), and p a sum of "localized charges" (Dirac delta 
functions). The charges are placed at the cone points' locations, thus corre- 
sponding to the irregular vertices of the mesh. The charge strength is equal 
to the cone excess angle (the difference between the cone angle and 2ir) which 
corresponds to the number of cells incident on the irregular vertex. For exam- 
ple, in the manifold shown in Fig. 2, (hi), the cone point has an excess angle 
of — 7r/2. The charge of the cone point, located at a point p, corresponding 
to such a excess angle is — 1 where 5^f> is the Dirac delta function in two 
dimensions. This is the charge of any singularity corresponding to an irregular 
vertex surrounded by three cells. Along with the boundary alignment condi- 
tions on 0, the problem of finding an appropriate manifold is thus reduced 
to an Inverse Poisson problem, of finding a charge distribution adhering to 
these conditions. The reduction gives exact, global relations that the function 
must obey. 

The inverse Poisson problem on planar domains is of interest in many fields of 
science and engineering (see references in section 7.3); the relevance of existing 
techniques to the present application remains to be examined. Another possi- 
ble application is to the problem of creating a surface mesh aligned with pre- 
defined directions, which has recently attracted much attention [4], [8], [13], [14]. 
The algorithm described in [13] creates conformal parametrizations and meshes 
approximately aligned with predefined directions. In that work, the local size 
and the local direction are treated as independent variables. In order to re- 
duce the number of singularities, a preprocessing step that modifies the cell 




Fig. 3. Boundary alignment. Boundary marked with heavy line. 
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size demand (the "curl correction" process) can be applied. In contrast, the 
theory developed in the present work uses conformality to a-priory link the 
cell-size and cell-direction fields, leaving just one field to work with, either the 
cell-size or the cell-direction. The applications of this theory to surface param- 
eterization and surface quadrangulation problems are an interesting subject 
for future work. 

Conformal parametrization of manifolds with conical singularities have been 
used in surface parametrization problems. For reviews and related work, see 
[9]- [13]. Mesh generation with boundaries and surface parametrization are 
different problems, because of the boundary alignment requirement in mesh 
generation. For example, the parametrization problem for planar domains is 
trivial: the coordinates of the plane form a good parametrization, but do not 
solve the mesh generation problem. 

The rest of the article is organized as follows: In section 2 some elements 
of differential geometry are shortly reviewed. In section 3 the cross-field and 
0-manifold are defined. Section 4 discusses the relation of the definitions in 
previous sections to mesh generation. Cone points are analyzed in section 5. 
The necessary and sufficient conditions for cross-field existence are developed 
in section 6. The set of conditions derived in sections 5,6 forms the core result 
of the work. In section 7 the theory is discussed through four case studies. Case 
Studies A,B are used to discuss the meaning of the conditions derived before. 
In Case Study C the possible structure of a mesh generation algorithm is 
discussed. An example quadrilateral mesh problem, though admittedly simple 
and artificial, demonstrates how, given an input, a manifold with cone points 
is constructed, adhering to the imposed conditions. This manifold allows the 
construction of meshes with irregular vertices; at any given point in the domain 
that is not an irregular vertex, at the limit of increasingly finer meshes, the 
cells' shapes tend to a square. Case study D gives an example of a curved 
surface meshing problem, and how it is solved. 



2 Parallel Transport and Geodesies 

In this section basic facts from differential geometry, required in subsequent 
sections, are shortly reviewed. More complete accounts can be found in any 
textbook on differential geometry, such as [15], [16]. 

For a surface D embedded in three-dimensional space, distances and angles 
on the surface can be defined by the embedding. In some coordinate system, 
denote by gij the metric given by the embedding. For example, on the plane 
with Cartesian coordinates, = 5ij, the Kronecker delta. Another metric, 
g^, is said to be conformally related to g^, if there exists a real function <p on 
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Fig. 4. Parallel transport of a vector G along a curve a. 
the surface such that 

9H = e 2 V (1) 

Given two vectors x = (x 1 , x 2 ),y = (y 1 , y 2 ) defined at some point on a surface 
or Riemannian manifold, the angle between them is given by 

cos# = (x'gijy* )/(yJ x k g kl x l ^y m g mn y n ) , (2) 

where summation over repeated indices (the Einstein convension) is assumed. 
The angle cos 6 as measured with the metric is 

cos 9 = {x% jy j )/(y/ x k ~g kl x l J y m g mn y n ) . (3) 

Substituting g^ = e 2 ^gij in Eq. (3 ), and comparing with Eq. (2) it is found 
that cos# = cos#, so the measurement of angles using the two conformally 
related metrics agrees. 

Given a vector field on the plane (with the Euclidean metric) the question of 
whether two vectors at two different points are parallel has a definite answer. 
This is not the case on curved surfaces, and more generally, for Riemannian 
manifolds. There, the definition of parallel vectors at different locations gener- 
ally depends on the path chosen between the points. For a curve a connecting 
points a and b, the parallel transport of a vector from a to b along a can be 
defined. If G is such a vector at a we denote its parallel translate to b along 
a by PT aJUb G, see Fig. 4. 

The geodesic curvature K g of a curve is the amount by which a curve "turns" . 
On the plane (with the Euclidean metric) turning is measured by the change 
of the angle of the tangent vector to the curve. On a surface (and, more 
generally, on a Riemannian manifold) the angle is defined relative to a vector 
that is parallel translated along the very same curve. For a curve a, denote 
the tangent vector at x by T a (x). Define 9 (x) — Z. (PT a a >x Tu (a) , T a (x)j. 
Then 

* = (4) 
where s is the length parameterization of a, see Fig. 5. 

The parallel transport depends on the metric, so a curve can have different 
geodesic curvatures under different metrics. Let K g , k g be the geodesic cur- 
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Fig. 5. Definition of n g . 

vatures of a curve a at some point x, with the metrics gij,(jij respectively, 
related as in Eq. (1). (Henceforth, all quantities relating to the metric (jij will 
be marked with a tilde.) Then related by 

k g = e"* (k 9 - d n 4>) , (5) 

where d n (f> = d(f>/dn is the derivative of along the normal vector N a , which 
is defined such that (T a , N Q ) form a right hand system. Eq. (5) is derived in 
appendix A. 

A geodesic is a generalization of a straight line. It is the curve whose geodesic 
curvature vanishes. The equation for a geodesic of the metric g^ is found by 
substituting k g = in Eq. (5): 

Kg = d n <p. (6) 



Two integral theorems are used throughout the paper. The first is the Gauss- 
Bonnet theorem, which relates the change in a vector undergoing parallel 
transport along a closed loop, to the total Gaussian curvature inside the 
loop. Let a(s) be a closed curve, a < s < b, a (a) = a (6), enclosing a re- 
gion R. The junction angles, 9i..9 N) measure the change in direction of the 
tangent at junction points. As a convention we assume that a is traversed 
in a counter-clockwise direction, let V be a vector at a (a). Then the angle 
I (V, PT aJU W) is equal to 

N 

I (V, PT a ^ a V) = 2tt - j>K g ds -J> = J J Kda (7) 
where K is the Gaussian curvature. 

The second integral theorem is Green's theorem, also known as the divergence 
theorem, or Gauss' theorem. Suppose a (s) is a curve enclosing a region R, 
transversed in a counter-clockwise manner, and a function defined on the 
surface. Then 

J J R V 2 cj)da = -j>d n <j>ds. (8) 
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The (unconventional) minus sign appears because the normal direction to 
a was defined such that (T a , N a ) form a right hand system, so N Q points 
inwards. For a surface the Laplacian operator V 2 denotes the Laplace-Beltrami 
operator [16]. 

If two metrics are conformally related as in Eq. (1 ), and the manifold with 
the metric is flat, i.e. 

K = 0, (9) 

a differential equation for can be derived. This can be done using the ex- 
pression for the curvature K in terms of It is derived here in a different 
way, using the integral theorems quoted above, as this technique is used again 
in subsequent sections. 

Suppose that a region R equipped with the metric is flat, i.e. Eq. (9) holds 
in R. Then according to the Gauss-Bonnet theorem, Eq. (7), 

N 

fk g ds = 2ir-^2ei. (10) 
The length element for the two metrics discussed are related by 



ds = \Jdx % gijdxi = \J dx l gijdxi = e^ds. (11) 
Substitute Eq. (5), (11) in Eq. (10) to get 

jk g ds = je~* (Kg - d n <p) e^ds = j>K g ds + J J W 2 <pda. (12) 

a a a 

Here Green's theorem, Eq. (8), was used. Subtracting Eq. (10) from (12), 

N 

= 2?r -J29i- fitgds- J J W 2 (j)da = j j (k - V 2 </>) da, (13) 
i=i a ^ 

where Gauss-Bonnet was used. Since this result is true for an arbitrary flat 
region R, the integrand K — V 2 must vanish, i.e., for any point in D where 
Eq. (9) holds: 

V 2 = K. (14) 
Eq. (14) is a differential equation for 0. It is a well-known result in conformal 
geometry, see e.g. [17], [18]. 



3 Cross-field and 0-manifold definitions 



The input to the mesh generation problem is assumed to be a surface D 
embedded in three dimensional Euclidean space, such that: 
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(i) Its boundary 3D is a union of a finite number of connected components 
' dD = UTj. 

(ii) Every component Yj is a piecewise C 2 closed curve. 

The points where a boundary component Yj is not differentiate will be called 
junction points. The set of all junction points will be denoted by J. 

The edges of the final mesh are to be laid along geodesic curves of a manifold 
that will be defined below, the -manifold. To create a high quality mesh, these 
curves should cross each other, and reach the boundary, at certain angles. In 
the case of a quadrilateral mesh, a cell with right inner angles is preferred for 
many applications. In the case of a triangular mesh, the preferred inner angle 
is 7r/3 radians. It is therefore natural to define the direction of edges at a point 
to within an addition of 71/ 2 radians in the quadrilateral case, and n/3 radians 
in the triangular case. This difference between the two cases leads to slightly 
different results; for clarity of presentation, the quadrilateral case is presented 
first, and the results for the triangular case are defered to Appendix B. 

The directions of edges at a point will be called a cross. The field of edge 
directions will be called a cross-field. 

Definition 1 (Cross) First, define an equivalence ~ of vectors in R 2 : for 2 
vectors vi, v 2 G R 2 ; Vi ~ v 2 if and only if Vi,V2 are parallel or perpendicular. 
A cross is an element 0/R 2 / ~. 

Thus, a given cross is a set of vectors, every pair of which are either perpen- 
dicular or parallel. 

If two vectors, defined at some point, undergo parallel transport along the same 
curve, the angle between two vectors is preserved: the angle before the parallel 
transport is equal to the angle after the parallel transport [15]. Therefore, the 
cross equivalence class structure is preserved, and the parallel transport of 
crosses is well-defined. 

Let P be a finite set of points in D = D U 3D. Denote the metric on D by g^ . 
Supoose a metric is defined on D\P. The following definition of a cross- 
field assures that the flow lines of the cross-field are geodesies, and that the 
crosses on the boundary are aligned with boundary. 

Definition 2 (Cross-field) A cross-field on a given manifold with metric 
is a mapping V : D\ (P U J) — > R 2 / ~ such that: 

(i) For points a,b G D\ (P U J), the parallel transport under the metric g^ 
of V (a) to b along a curve a is independent of a, and is equal to V (b): 
PT a ^ b V(a) = V(b). 

(ii) For a G Tj\(J U P), the tangent belongs to the cross there: Tjv (a) G 
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Fig. 6. Illustration of a cross-field, and its cross-field geodesies. The cross-field is 
represented by the crosses. Thin lines represent the cross-field geodesies. Heavy 
lines represent the boundary. Crosses at points g, e are shown. Points c, k, b are in 
P, hence no crosses are defined at these points, h is a junction point: h G J. 

V(a). 

The flow-lines curves of a cross- field are geodesies of the metric g^. That is, a 
geodesic aligned with the cross-field at one point (i.e., whose tangent belongs 
to the cross at that point), is aligned with the cross- field everywhere else on 
the curve. This is because both the cross-field and the tangent to a geodesic are 
parallel-translated along the geodesic, see definition 2,(i), and the discussion 
preceding Eq. (6). Geodesic curves aligned with the cross- field will be called 
cross-field geodesies. Cross-field geodesies and their relation to the cross-field 
are illustrated in Fig. 6. 

As before, let P be a finite set of points in D. We now define the (p-manifold. 

Definition 3 (0-manifold) A 0-manifold is a Riemannian manifold defined 
on D\P , with metric (jij such that: 

(i) gij is conformally related to g^, i.e. g^ = e^gij where is a real function 
on D\P. 

(ii) The <p-manifold is locally flat, i.e. its Gaussian curvature tensor K = 
for all points in D\P. 

(Hi) For every boundary point a G Tj, a J U P, the limits lim r ^ a (p (r) and 

lim r ^ a V0 (r) exist and are continuous along Tj at a. 
(iv) A cross-field with the metric g^ exists. 

The points in P will be called cone points. This name is justified in section 
5, see also section 1. For now, the points in P are just points where is 
undefined. 

In the following sections, the requirement that a 0-manifold exists is translated 
into conditions on the function 0. 
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Fig. 7. A geodesic "rectangle" . 
4 Relation to mesh generation 

The cross-field formulates the demand that mesh cells have certain inner an- 
gles. In order to have square-like shapes, the cells should further have edges of 
similar lengths. This is where the conformal metric property of the 0-manifold 
comes into play. 

Integrating over Eq. (11), the length of a curve a on a 0-manifold is given by 

s(a)= [ ds= [ e^ds. (15) 

The area of a region of a manifold is given by [15], [16] 

A(R) = J J R Jgdx 1 dx 2 = J J R e 2 ^^dx l dx 2 , (16) 

where g = det (gij) , g = det The tilde denotes, as before, a quantity 
with respect to 0-manifold metric gij. 

The following claim is a consequence of the isometry of a flat manifold to 
the Euclidean plane [16]. It is the "manifold version" of the properties of a 
rectangle. 

Claim 1 Let 71, 72, 73, 74 be 4 geodesic segments in a flat region of a manifold, 
organized as in Fig. 7. Suppose that the inner angles at the vertices A, B, C are 
right angles, and s (74) is the manifold-length of the i-th side of the "rectangle" . 
Then the inner angle at D, Od, is a right angle, and s (71) = s (72) , s (73) = 
s (74). The area of the "rectangle" is s (71) • s (73). 

On an Euclidean plane, the properties of a rectangle allow one to lay a grid 
on a region of the plane. A grid can be regarded as two families of mutually 
perpendicular straight lines, with equal spacing between lines of each family. 
The same can be done for a 0-manifold using the "rectangle" properties stated 
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Fig. 8. A grid on an Euclidean plane (left), and on a flat manifold (right). 

e _< ^ large 




Fig. 9. Geometric interpretation of the formula n g = d n (j): manifold-edges curve 
towards smaller cells. Cell size is proportional to '. 

in claim 1, see Fig. 8. The grid divides the space into square-like regions, 
each bounded by four geodesic segments of manifold length As, that will 
be called manifold- edges. The regions enclosed by the manifold-edges will be 
called manifold- cells. 

As is a single, fixed number for the mesh. The overall size on the surface or 
plane of the manifold-cells can be controlled by adding a constant to 0, so 
the freedom in choosing As is redundant, and we set As = 1. According to 
Claim 1 such a cell has unit -manifold area. According to Eq. (15), for small 
enough manifold cells (large enough 0) , a manifold edge of length As = 1 has 
length As (j e ) ~ e~^, and e~^ is interpreted as the local edge length, or local 
cell size. 

A geometrical interpretation of the relation k 9 = <9„0, Eq. (6), can now be 
given. For a mesh with approximately square cells, the curving of manifold- 
edges is related to the changes in cell size in the perpendicular direction, see 
Fig. 9. This is quantified in Eq. (6): k 9 is the curvature of the lines defining 
the edges, and <9 n is the change in 0, which is related to local cell size by e~^. 
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Fig. 10. Parallel transport around a cone point. 



5 Cone-points 

The existence of a cross-field restricts the possible behavior of 4> in the vicinity 
of cone points. Let p £ P be a cone point, and a £ D\P a point that is not 
a cone point. Let a be a simple closed curve starting from a that encloses p, 
and only p of P, see Fig. 10. The junction angles of a are denoted by ai-.a^. 
According to cross-field definition, the cross V (a) is parallel translated to V (a) 
along a, that is, a vector y £ V (a) is parallel translated along a, to y' £ V (a). 
By the defintion of a cross, y, y' are either parallel or perpendicular, and the 
angle Z (y,y') is kir/2, k £ Z. Using Eq. (5), (7): 



where 5 is the region enclosed by a. 

On a 0-manifold a manifold-cell has unit area, so in order to have a finite 
number of cells, the the area of the 0-manifold must be finite. This will further 
restrict the type of singularity allowed at a point p, as is now shown. 

For a cone point p on the surface let U be a neighborhood of p in which 
there exist isothermal coordinates (x 1 ,^ 2 ). In such coordinates, which can 
always be found localljdl the metric takes the form: = F~ 1 5ij, where 
^(x 1 ,^ 2 ) is a real function of (x 1 ,^ 2 ). Note that on the plane, the standard 
Cartesian coordinates satisfy g^ = 8^, and hence are isothermal. In isothermal 



2 The original proof, due to Gauss, requires that the surface is analytic [16]. There 
are proofs with weaker assumptions, but these distinctions are immaterial for the 
present purposes. 




(17) 
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coordinates the Laplace-Beltrami operator can be written a 

V 2 = W*\* 2 )(^ + ^V (18) 
V ' \d(x^) 2 d(x 2 ) 2 ) V ; 

Eq. 18 can be seen as a generalization of the usual planar Laplacian, given 
by F (x l , x 2 ) = 1. Let Br be a disk of radius R in (x 1 , x 2 ), i.e. the region for 
which (a; 1 ) 2 + (x 2 ) 2 < R 2 . In isothermal coordinates Eq. (14) reads 

/ d 2 (h d 



\d (x 1 ) d (x 

in the punctured disk Br\ {p}. In polar isothermal coordinates (r,ip), (that 
is, the polar coordinates corresponding to (x 1 ,^ 2 )), the general form of this 
solution can be written as (see e.g. [20]): 

Q oo 

(r, -0) = 4>p + — In r + V b n r n sin (nip + c n ) , (19) 

27r n=l 

where <pp is a solution to the Poisson equation V 2 0p = K/F in Br (including 
p), and Q, {& n }^Li , { c n}^Li are all real numbers. 

Let a be a curve tracing the circle (x 1 ) 2 + (x 2 ) 2 = R 2 , in the counter-clockwise 
direction. For given by Eq. (19), the flux of V0 through a is 



- y -^ds = J J B K/F^gdx l dx 2 + Q, (20) 

where Green's theorem, Eq. (8) was used to convert the first term to a surface 
integral. In the limit R — > Eq. (20) becomes: § a ^ds — > —Q, and Eq. (17) 
with S — > becomes: f Q f^rfs — > /cV/2 for some k' G Z. Comparing the two 
limiting values for § a ^ds we find 

Q = k\ (21) 

for some fc 6 Z. 

The quantity Q will be called the charge of the cone-point. Eq. (21) states 
that the charges must be multiples of n/2. 

In order to have a finite number of unit-area manifold cells, the manifold-area 
of a neighbourhood of p must be finite. The manifold area of the disc Br of 



3 This can be seen by substituting the form of the metric tensor in isothermal 
coordinates g^ = F~ 1 5n : into the definition of the Laplace-Beltrami operator: 
V 2 (/> = g ik (j)n-,k, see e.g. [16]. 
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(i) 



(ii) 



(iii) 






Fig. 11. Cross-field geodesies around a singularity. The geodesies are placed at unit 
manifold-distances apart, (i) A k = — 1 singularity, (ii) A k = 1 singularity, (iii) A 
k = —4 singularity, with an infinite number of cells. 

radius R around p is obtained by substituting Eq. (19) into Eq. (16): 



A(B 



R 



e 2 ^ '^fgdx l dx 



U~2 



B f 



r R ( k 00 \ 

/ exp 20p H — lnr + 2 2J ™ sin (n^ + c n ) y/g2Tirdr 
Jo \ 2 n= i ) 



2tt / r 
Jo 



k/2+1 



exp 2 J] b n r~ n sin (n^ + c n ) e 2<t>p Jgd 



[22) 



n=l 



Since a solution of V 2 0p = fT is bounded in a compact region, e 2 ^ p does not 
effect the convergence of the integral, nor does ^/g, which is bounded away 
from zero. To converge, it is required that b n = for every n, and that k > —4. 
This restricts <b to the form 



k 

(r) = cf)p + — lnr. 



(23) 



with k > —4. This is a solution of the Poisson equation for one point-source 

(24) 



V 2 <P = K + ^-5i 2 \ 



in a neighborhood of r = 0. For many cone points Eq. (24) becomes: 

Condition 1: 

(p obeys the equation 

V 2 <P = K+l £ kiS® (25) 

Z i=l..N 

the Poisson equation with point sources (delta functions) 5^ , with fc, E Z, 
h > -4. 



For a planar domain, Condition 1 can be rewritten as 
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Condition 1 (planar domain): 



For r G D\P, (f) (r) can be written as 



<p{r) = <p L + ] - ki\n(\r -pi\) , 



i=l..N 



with ki > —4, P = {pi} i=1 N , and V 2 0l = on D. 

Fig. 11 shows selected cross-field geodesies for a planar domain around a singu- 
larity (with (ftp = 0), spaced at manifold-distances of As = 1 from each other, 
for different singularity strengths. In practice, good candidates for meshes will 
only have singularities of charge k > —1, due to the inner- angles of cells 
incident on the singularity, see also Remark 7 in section 7.3 below. 

To summarize this section, it has been shown that the existence of a grid of 
geodesies that follow a cross-field and create a finite number of cells, restrict 
the form that the function can take. This is formulated in Condition 1, 
stating that <fi must obey the Poisson equation, with delta function charges 
corresponding to cone points. 



6 Boundary Alignment Conditions 

This section examines the conditions for boundary alignment of a cross-field, 
see Definition 2,(ii). Three conditions will shown to be necessary. Sufficiency 
of the conditions is discussed in section 6.1. 

We start by developing an equation that will be used in the derivation of 
the conditions below. Let a±, 02 be two points on two boundary curves r\, r 2 , 
respectively. Note that ri,r 2 may be the same boundary curve: Ti = T 2 . Let 
a be some curve from a\ to a 2 , see Fig. 12. For % — 1, 2 denote by 9 ai = 
Z (Tr, (aj) , T a (aj)). By definition of a cross-field, the cross in a\ must be 
parallel-translated along a to the cross in a 2 . The cross in a\ contains Tr x (ai), 
and the cross in a 2 contains Tr 2 (a 2 ) so 




^ (PT ai ^ a T Tl (a,) , T r2 (a 2 )) = rj + 9, 



(26) 



where 77 is defined as 



V = I (PT ai ^ a T ri (aO,T Q (a 2 )) 



= ^ ( PT ai ^a 2 T ri ( ai ) , PT ai ^ a T a ( fll )) + I {PT ai ^ a T a ( 0l ) , T a (a 2 )) 
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Fig. 12. Condition 4. 

a 




Fig. 13. Condition 2. 

The last equality is due to the preservation of angle in parallel transport, 
together with Eq. (4). Eq. (5) and (11) can now be used with Eq. (27) 



Eq. (28) is used below to develop the next three conditions. 

Let (3 be a segment of a boundary curve Tj, between points a,b G Tj that are 
not junction points: a,b ^ ( J U P), and suppose (3 does not contain junction 
points. Choose some smooth curve a from b to a, whose tangents at a, b coin- 
cide with the tangents to f3, see Fig. 13, i.e. 9 a2 = 9 ai = 0, and such that the 
region enclosed by a, f3 does not include any cone points. According to Eq. 




Substituting this into Eq. (26) we find the relation 




(28) 



(28), 




(29) 



According to Gauss-Bonnet, Eq. (7), 




(30) 



where S is the area enclosed by a, (3. Green's theorem gives 




(31) 
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Substituting Eq. (30), (31) and recalling that V 2 = K in D\P (see Eq. (14)), 
Eq. (29) now reads 

nir/2 = —2n + / K g ds — / d n (j)ds. (32) 

J/3 J/3 

When a — > b, the integrals both tend to zero, and the equality is possible only 
if nn/2 = —2ir. Then 

/ (k 9 - d n (p) ds = , (33) 

J/3 

and with a — > b it follows that 
Condition 2: 

For a point a G Tj\ (J U P) boundary curvature K g , <p must satisfy 

d n (f) = Kg. 

The equation in Condition 2 is the same as Eq. (6). This is not incidental: 
Condition 2 causes cross-field geodesies, that are close and parallel to the 
boundary, to follow the shape of the boundary, as shown in Fig. 3. 

For a point c G (J U P), that is, a junction or cone point of the boundary, let 
a, b G Tj\ (J U P) be two points on both sides of the junction point c, and (3 
the boundary segment from b to a, see Fig. 14. Let a be a curve from a to b. 
The Gauss-Bonnet theorem, Eq. (7), for the curve [a,/3] reads 

/ Kgds = 2ir — Kgds — 9 a + 9 h — (it — 9 in ) — Kda, (34) 

J a J fi J J S 

where 9 in is inner angle at c, and 9 a ,9 b are defined as in Eq. (28). Eq. (28) 
reads 



9 b ~9 a = / (n g - d n <p)ds + n-. (35) 

J a Z 

Adding Eq. (34) and Eq. (35) and rearranging: 

J d n (j)ds = n — + 9 in — J Kgds — J J Kda. 

When a — > c,b — > c then fp K g ds — > 0, / J* s Kda — > 0, leading to: 
Condition 3: 

For a cttry e a as described above: 

f 71 

/ <9 n 0ds = n- + 9 in , (36) 
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Fig. 14. Condition 3. 




Fig. 15. Condition 4 for a planar domain. 

with n G Z. This requires <9„0 to have a singularity at c. On the plane, for 
example, <p has a singularity of type |^ ~ 1/r at distance r from c. 

The final condition to be imposed on is a relation between different boundary 
components, and accordingly, its formulation is not local. It is just a restatment 
of Eq. (28). 

Condition 4: 



For two boundary components ri,r 2; and a curve a connecting ai G Ti, to 
d2 G T 2 it is necessary that 

rai ra.2 7f 

/ d n (f)ds = 6 a2 -6 ai + / K g ds + n-, (37) 
for an integer n. 9 ai ,9 a2 are defined as in Eq. (28). 

Condition 4 states that the total flux through a curve connecting the two 
boundary components can only belong to a certain, discrete, set of values. 

Condition 4 has can be put in a simpler form when the domain is planar. Let 
Ot be the angle between the x-axis and T Fl (ai), see Fig. 15. Then, because 
lai K gd s is equal to the change from a x to a 2 in the angle between the tangent 
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to a and the a>axis we have 




/ 

J a 



Kgds — 9 2 — Oi- 



Condition 4 for the planar case then reads 
Condition 4 (planar domain): 

For two boundary components T 1 ,T 2 , and a curve a connecting a x G r 1; to 
a 2 G T 2 it is necessary that 



for an integer n. 

6. 1 Sufficiency of the conditions 

In the previous sections Conditions 1-4 were shown to be necessary for the 
existance of a 0-manifold. This section considers the question: when are Con- 
ditions 1-4 sufficient? It turns out that the answer to this question depends on 
the genus of the surface. Intuitively, the genus of a surface is the number of 
handles a surface has [22]. For example a sphere or a disk are genus-0 surfaces, 
and a torus is a genus-1 surface. A surface of any genus can have any number 
of boundaries: for example, any planar domain, with an arbitrary number of 
boundaries, is a genus-0 surface. 

The following theorem states that for genus zero surfaces, Conditions 1-4 are 
sufficient. The theorem also shows how to construct the cross-field given the 
function 0. For higher genus surfaces, an additional condition is required, 
constraining the parallel transport along curves "passing through" the handles 
of higher genus surfaces. A detailed disscusion is out of the scope of the present 
work, but a brief account is given in Appendix C. 

Theorem 4 Suppose D is a surface of genus-0. Let ao G Tj be some boundary 
point such that a ^ J U P. Assume <fi satisfies Conditions (1),(2),(3), and 
satisfies Condition (4) with curves between a and points G r i; for every 
i jo- Let V (ao) be the (unique) cross such that Tr JQ (a G ) G V (ao)- For every 
b G D\P, let a be some curve from a to b, and define V (6) to be the translate 
of V (ao) along a. Then V (6) is independent of a, and V is a cross-field. 

The proof is given in Appendix C. 

Remark 5 If the surface is closed, i.e. has no boundaries, the cross at some 
point a on the surface must be fixed for the cross-field to be unique. 
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7 Discussion: Finding a 0-manifold 

According to Theorem 4 in section 6.1, if a function is found, satisfying 
Conditions 1-4, a 0-manifold exists, and its cross-field can be calculated. 

In addition to Conditions 1-4, a meshing problem can include other require- 
ments, such as cell size requirements on the boundaries, or inside the mesh. 
Clearly, finding an appropriate 0-manifold then depends on these conditions 
as well. In what follows the problem of finding a 0-manifold under different 
requirements will be discussed. 

The next three sections focus on the planar case, and give concrete examples 
of the theory presented above. The fourth section gives an example of a curved 
surface meshing, that is solved analytically. 

After finding a valid 0-manifold, the final stage of a mesh generation process 
involves finding a discrete partition into well-shaped manifold cells, see below. 
A comprehensive discussion of this step will not be given here, but some 
comments on this process will be made as part of the case studies. 

1.1 Case Study A: No singularities are needed 

Sections 7.1,7.2,7.3 deal with meshing of planar domains. On the plane, the 
geodesic curvature k 9 reduces to the curvature of a planar curve, which will 
be denoted by k. 

We start with a simple planar case. Suppose that the boundary has only one 
connectivity element (loop) F 1 . Furthermore, suppose that all junction angles 
6ji (equal to tt — 6 in , 6 in the inner angle) are multiples of 7r/2, and that the 
their sum is J2i ®.n — 27r. In such a case, Condition 4 is empty (since there 
is only one boundary loop), and Condition 3 can be satisfied with = — 1 , 
and without any additional singularities at the junction points. Condition 2 
are Neumann boundary conditions on <fi. T\ is a simple loop, and the sum of 
junction angles is 27r, thus the total is flux of through the boundary is 



Therefore given the boundary conditions a solution to the Laplace equation 
(that is, with no singularities) exists, and is unique up to an additive constant 



A simple example of such a case, that can be solved analytically, is a section 
of an annulus between two radial lines, see Fig. 16. The boundary conditions 




[19]. 
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Fig. 17. C-frame geodesies for the example in Fig. 16. (i) equally-spaced geodesies, 
(ii) geodescis forming well-shaped manifold cells. 

on the sides of the boundary, dictated by the shape of the boundary are given, 
according to Condition 2, by (note that k can be negative, see Eq. (4)): 



d<j> 
dn 


i 


= Ki = 


1 

= % 


(38) 


d<f> 
dn 


2 


= K 2 = 


1 




d(j> 




def) 


= 0. 

4 




dn 


3 


dn 





The solution to the Laplace equation V 2 = with boundary conditions given 
in Eq. (38) is 

(j) = -\nr + C. (39) 

Here, r is the distance from the center of the ring, and C is a constant. Note 
that no specification of the local cell size was given in the input, and, according 
to section 4, the resulting local cell-size function is part of the solution. It is 
given by exp (—0) = Cr, with C = — InC". 

Fig. 17 shows geodesies starting from the boundaries. To create the figure, as 
well as Fig. (18), Fig. (19) and Fig. (20 ),(hi), the Poisson equation was solved 
numerically on a triangular mesh inside the domain (even though an analytical 
solution is known in the case shown in Fig. (17)). The geodesies were calculated 
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Fig. 18. A more elaborate exmaple than that of Eq. (38). (i) Flow-lines of V</>. (ii) 
Equally-spaced geodesies. 

by solving the geodesic equation, Eq. (6), with initial tangent perpendicular to 
the boundary. The integration constant C was chosen such that 10 manifold- 
edges will fit on the radial boundaries, i.e. that /j? 2 e^dr = 10. This fixes the 
solution completely, and gives a non-integer manifold-length along the arcs. 
Fig. 17, (i) shows geodesies spaced at unit manifold- distances of each other. 
Whereas the radial direction fits exactly 10 manifold-edges of equal length, 
the manifold- distance from the left-most geodesic to the boundary is less then 
1, resulting in manifold cells with high aspect ratio in the left-most row of cells. 
In Fig. 17, (ii), two different spacings are used, so as to allow equal spacing 
in both the radial and tangential directions, with an aspect ratio as close to 
one as possible. Fig. 18 shows an example with a more elaborate boundary 
adhering to the restrictions stated in the beginning of section 7.1. Again, the 
equally-spaced geodesies of Fig. 18, (ii) do not form well-shaped (or even valid) 
cello. A valid discretization can be created, e.g. by using geodesies emanating 
from the junction points to decompose the domain. 

Remark 6 Note that in general, local edge directions are not aligned with the 
flow lines of V0 ; as can be seen e.g in Fig. 18. 

We now lift the restriction that all junction angles are multiples of tt/2. If 
the angle is not n/2, Condition 3 requires that <fi have a radial-singularity at 
the junction point. Denote the junction inner-angle by 6 in . Suppose that 
contains a singularity caused by placing a charge at the junction point, i.e. 
<fi = Q- In r. (If two or more boundary segments are incident on the same point, 
other functions may be required.) Let a r be the curve formed by traversing an 
arc of the circle at a distance r from the junction point in a counter-clockwise 
direction, as in Fig. 14. Then according to Condition 3 (Eq. (36)) 

—n + Ui n — / — as — — - — U in r — 

2 Ja r on 2nr 27T 



4 The geodesies shown near the right and bottom of the intrusion in Fig. 18 are 
not parallel to the boundary. This is due to the change in cell size. Other geodesies, 
closer to the boundary, follow the shape of the boundary more closely. 
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for ml G Z, so 

Q = 2vr (n^? - l) (40) 

for n = —n'. n must be positive, since otherwise Q > 2tt, causing the manifold- 
area to diverge at the singularity, by the same argument as presented in sec- 
tion 5. Apart from this restriction the number n is not fixed, and affects the 
manifold- angle at the singularity, i.e. the number of manifold-cells incident on 
the junction. 

(i) (ii) 




3 / 



/ 



Fig. 19. A domain with a singularity at a junction point. 

The example in Fig. 19 shows a domain enclosed by one boundary element. 
In one junction, the inner angle is 7r/3, which, according to Eq. (40), requires 
a singularity of charge Q = n, for k — 1 placed at the junction. The solution 
was calculated by decomposing <fi into 2 contributions: = C + 0l. C is 
the charge potential, C = ^-lnr. (f>i was computed by numerically solving 
the Laplace equation on a triangular mesh with boundary conditions = 

d<t> d<t> c 

dn dn 

The examples that were presented in this case study could have been obtained 
by a conformal mapping from a "logical" domain. This is, of course, not true 
when cone points are present inside the domainal. 

Note, however, that unlike many mapping techniques, even if such a "logical" 
domain can be defined, its shape is not fixed in advance, and is part of the 
solution. This is even more pronounced in problems involving cone points, see 
below. Conformal mappings that are also boundary aligned are quite limited 
in the scope of problems they can mesh, and sometimes yield large differences 
in cell size (as in the example shown in Fig. 18). That is why in mapping 

5 Even without cone-points inside the domain, this is not always possible, since the 
mapping from D to the "logical" domain is not, in general, one-to-one. In manifold- 
theory terminology, even if the manifold is flat and simply connected, it is not 
necessarily covered by a single geodesic coordinate patch of the conformal metric. 
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(i) 




o 




Fig. 20. The significance of Condition 4. (i) The input boundary, (ii) Ignoring Con- 
dition 4, <j) = const is a possible manifold. The resulting cross field does not con- 
form with all boundaries, (iii) A manifold with two singularities. Selected cross field 
geodesies are drawn. 



techniques the conformal restriction is lifted, see e.g. Ref. 2. In the present 
work, the conformal condition (in its manifold formulation) is retained, and 
instead cone points are allowed into the manifold. 



7.2 Case Study B: Two boundaries, no boundary cell-size demand 



The purpose of the example in this section is to demonstrate the meaning and 
relevance of Condition 4. Unlike the other boundary alignment conditions, the 
formulation of Condition 4, stating the required relations between different 
boundary elements, is not local. We now show that the relative placement of 
different boundary elements can force the introduction of cone points in order 
to obtain a valid cross-field. 

Fig. 20, (i) shows a domain bounded by two boundary elements, an "outer 
loop" Ti, and an "inner loop" r 2 . Assume that there are no cell-size demands. 
We start by ignoring Condition 4, and trying to proceed as in Case A, that 
is, searching for a solution without any cone points. The curves composing 
r!,r 2 are straight lines, so the Neumann boundary conditions read |^ = 0, 
and the solution to the Laplace equation is trivial: <fi = const. Condition 3 is 
also fulfilled with k — 1 at all junctions. However, this solution does not give 
a valid cross-field. A cross-field geodesic running from Fi to r 2 will not reach 
T 2 at a right angle, see Fig. 20, (ii). Thus, we cannot do without Condition 
(4), and since the only solution without cone points with boundary conditions 
|^ = is 4> = const, we learn that a solution without cone points does not 
exist. One possible solution with cone points is shown in Fig. 20, (iii). This 
solution, obtained using symmetry arguments, contains two cone points with 
opposite signs. Selected cross-field geodesies are shown in Fig. 20, (iii). One of 
them runs from one boundary to the other. The rest are geodesies that are 
incident on the cone points. 
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Remark 7 In Fig. 20, (in), cross-field geodesies reaching the cone points are 
shown. Three cross-field geodesies reach the left cone point, which has a charge 
of k = —1, and five reach the right cone point, which has a charge of k = 1. 
Moreover, it seems they reach the cone points at equally distributed angles. 
These observations are true in general. It can be proved that there are always 
exactly k + 4 angles from which cross-field geodesies are incident upon a sin- 
gularity, and these angles are equally distributed around the cone point. Such 
geodesies will be called star-geodesics. Note that this affects the angles of the 
mesh-cells created: for small cells with a cone point on one vertex, that vertex's 
inner angle will be approximately 7^7. 



7.3 Case Study C: boundary cell-size demand; finding cone points' locations 



In Case Studies A,B above no constraint on the size of boundary edges was 
given. Yet the boundary edge-length is often specified in meshing problems. 
Suppose we are given a function F stating the required cell size (that is, cell 
edge length) at each point on the boundary. The local cell size is e~^, thus 
F = e~^, giving a Dirichlet boundary condition on 0: 

0| r = -lnF. (41) 

Condition 1 states that obeys a Poisson equation with point charges playing 
the role of cone points. The problem of finding a suitable is reduced to 
finding a charge distribution (number of charges, their locations and charge- 
strengths), such that will fulfill Conditions 2,3,4, together with Eq. (41). 
This is an Inverse Poisson (IP) problem . As opposed to a Direct Poisson 
problem, where the charge-distribution p in V 2 = p is known, and one is 
asked to find 0, in an IP problem, certain information on is given, and the 
charge distribution p is to be found. 

IP problems have important applications in various areas of science and engi- 
neering [23]- [27]. By its nature, the IP problem is ill-posed, and the solution 
may not be unique, and may be sensitive to small changes of the input, such as 
small changes in boundary conditions. In problems of this type any prior infor- 
mation on the charge distribution can play an important role in the solution 
of the problem. 

The problem of finding can be broken into the following steps: 



(i) Given the boundaries Tj of the domain, and the cell-size requirement F 
on the boundary, calculate the Neumann boundary conditions = k 

J 1 J on Q£) 

(Condition (2)), and Dirichlet boundary condition 4>\ dD = — In (F) (Eq. 
(41)). 
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(ii) Impose boundary condition (3), e.g. by placing charges at the junction 
points (see Eq. (40)). 

(iii) Solve the IP problem: Find the (finite) number, location and strength of 
charges such that Neumann and Dirichlet boundary conditions calculated 
in (i), and Condition 4 hold approximately. According to Condition 1, the 
charges should be of strength fcj7r/2, with ki > —4. The charges can be 
placed: 

(1) Inside D (forming the set P). 

(2) At junction points (which amounts to changing k in Eq. (40)). 

(3) On the rest of the boundary (forming the set PddD), where according 
to Eq. (40) with in = ir, charges of strength kiii can be placed. 

(iv) Once the charges are placed, <fi is found by solving the standard (Direct) 
Poisson problem. 

Remark 8 (i) Note that since this is an inverse problem, i.e. the charge 

distribution is not fixed, both Neumann and Dirichlet boundary conditions 

can be imposed together, 
(ii) Though the limit on ki due to Condition (1) is ki > —4, charges should be 

of charge ki > —1 in order to have convex cells, and preferably with ki < 2. 

This is due to inner angles of cells incident on the singularity, see Remark 

7 in section 7.2. 



no 




Fig. 21. The steps of a possible mesh-generation process, (i) The input boundary 
(thick line). Locations of singularities of 4>q that was used to create the boundary 
are marked. (Square: k = 1, Circle: k = —1). (ii) A solution recovered using an IP 
algorithm (see text). Reconstructed singularities marked as in (i). Equally-spaced 
geodesies are shown (thin lines), (iii) Geodesies of the reconstructed solution, at 
approximately equal spacings. 

The steps outlined above are illustrated in Fig. 21 . For the purpose of this 
example, an input to the algorithm described above was created artificially by 
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joining four geodesies of a 0-field chosen in advance, at right angles to each 
other. The 0-field chosen for creating the boundary is a sum of fields from 
four charges: O ( r ) = Y,f=i l n l r — r i|> with the charges' strengths and 
locations. Two charges are located inside the boundary, and two outside, see 
Fig. 21 ,(i). This defines the shape of the boundary. The cell size requirement 
imposed on the boundary was F = e~^°, plotted in Fig. 22. The boundary 
shape and the cell size demand are the input of the problem. 




upper / lower boundary 

right boundary 

left boundary 



'0 1 2 3 4 5 6 

length along curve 

Fig. 22. The cell-size input requirement along the four sides of the boundary shown 
in Fig. 21, (i). The upper / lower boundary curves are traced from left to right. 

Using this input, the steps of the algorithm outlined above were followed. 
First (step (i)), Neumann and Dirichlet boundary conditions were calculated 
from the input. Step (ii) was fulfilled automatically without adding additional 
charges, because all junction have right inner angles. In step (iii) an algorithm 
solving the IP problem [28] was invoked. The algorithm exactly reconstructed 
the location and charge of the two charges inside the domain. (It is important 
to note that the IP algorithm of [28] could be readily applied due to the 
artificial nature of the example: the input was constructed such that these 
two charges will reconstruct the input conditions exactly. This may not be the 
case in other cases, and may require other IP solution methods.) The 0-field 
was then constructed by solving the direct, standard Poisson problem, with 
given charges. In this example, since the location of the charges inside the 
domain where recovered exactly, the (f) recovered was exactly <f)Q. 

Fig. 21, (ii) shows cross- field geodesies at equal manifold-distances As = 1. 
As was discussed in section 4, the function <p is defined up to an additive 
constant, that controls the overall cell-size. This constant was chosen such 
that the manifold-distance As between two geodesies incident on the two 
charges be equal to one, see Fig. 21, (ii). In Fig. 21 ,(iii) cross-field geodesies 
spaced at approximately equal manifold-distances are shown, such that no 
high aspect-ratio manifold cells, as those in Fig. 21, (ii), exist. Note that the 
star-geodesics divide the domain into sub-domains without charges, that can 
be meshed more easily. This hints at possible ways of using the 0-manifold for 
creating meshes. 
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Other mesh requirements can be formulated. Examples include cell-size within 
the domain (as in adaptive meshing), and cell direction. We comment briefly 
on these subjects. In the case of adaptive meshing, a requirement for cell size 
F as a function of location inside the domain is given. This is translated to a 
requirement on by = — In (F). This fixes completely, and therefore may 
not be fulfilled exactly along with other conditions, such as having V 2 = 
(in the planar case) almost everywhere. However, weaker constraints may be 
given, such as specifying F along curves within the domain. When specifing 
cell direction, the direction of the cross-field is given. This translated to a 
constraint similar to Condition 4, on the flux $ through some arbitrary curve. 
As is the case with adaptive meshing, a requirement for cell direction every- 
where inside the domain is too restrictive, but more limited demands may be 
applied, such as an approximate alignment. 



7.5 A curved surface example 



For many surface meshing applications it is desired that the cells be aligned 
with certain prescribed directions, such as the principle curvature directions 
of the surface. As an example of a meshing problem of this type, the class 
of surfaces known as surfaces of revolution is analyzed analytically in this 
section, and it is shown how a cross-field aligned with the principle directions 
is constructed in this case. 

A surface of revolution is defined by a curve in the (r, z) plane, that is revolved 
around the z-axis, see Fig. 23 . The direction of rotation will be called the 
^-direction. The line traced by the curve at a given 9 is known as a meridian. 
The circles at given (r,z) values are known as circles of revolution [15]. The 
principle axis directions at a point are the directions of the circle of revolution 
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and meridian passing through that point. A mesh aligned with these directions 
has iV cells around every circle of revolution at a given (r, z) . The cell size at 
a given point should therefore be N/ (2%r). But the cell size is equal to e"^, 
so 4> is expected to be 

# = -h(£). (42) 

We now show that this is indeed the case. By the alignment requirement, 
edges must be laid along circles of revolution and meridians. Therefore, these 
must be geodesies of the manifold with metric g^. Let T Q (s) be the tangent 
to the curve a (s) tracing a circle of revolution at (r, z), see Fig. 23. Then the 
geodesic curvature K g (s) is equal to 

dT a 1 

K„ = — • 1 a = - COS t) 

ds r 

where the angle 6 between a and the meridian, (3, is equal to cos 

so 

1 

Kg — 



V(^) 2+i 



For a geodesic of g^-, = — d(f)/dn. The normal is directed along the 

meridian, so can is found by integrating along the meridian (3. Noting that 

ds = \j Jj^j + Idr along (3, we have 



^ = / 7T ds = I ; ^ =~ I ~dr = - \n—. 



As expected in Eq. (42). The constant C is determined once for the surface, 
e.g. by the number of cells on a given circle of revolution. By comparing with 
Eq. (42), C is interpreted as C = N/2n. 

If the surface of revolution reaches the z-axis at some point p, and is not 
parallel to the z-axis there, the manifold will have a cone points of charge 
— 4| at p. To see this, we again use the fact that for a circle of revolution 
= k g = Kg — d(p/dn. Integrating over a revolution circle surrounding p: 

°=J( K °- 9 i) ds = 2 *-J 'I Kda+ J 7 v2 ^ tt 

= 21+/ \n-Sfda. 

Eq. (7), (25) were used. Thus, by / / 5^ da = 1 we have n = —4. We note that 
in the vicinity of a singularity of charge —4^ the number of cells diverges, 
see section 5. If this is undesirable, a mesh that isn't exactly aligned with the 
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principal directions may be constructed, e.g. by replacing the single charge 
with a few cone points with smaller charges. 



8 Conclusions 

In this work a continuum description of unstructured meshes was proposed. The 
structure aims at describing meshes that, away from irregular vertices, have 
well-shaped cells. In the limit of an increasingly finer mesh, the cell's shape 
approaches the shape of a square (in the case of a quadrilateral mesh) or equi- 
lateral triangle (for triangle meshes). Accordingly, in the continum limit such 
meshes can be described by just two local properties: the local cell size, and 
the local directions of edges, formalized by the notion of a cross. 

The connection between cell size and cell direction is established by defining 
a Riemmanian manifold, the 0-manifold. The geodesies of the manifold trace 
the edges of the mesh, as is formalized via the definition of a cross-field. This 
analysis allows the focus to turn to the irregular vertices, represented in the 
continuum structure by cone points. 

The demand that the mesh conform to the boundary, and have a finite number 
of cells, produces conditions on the function <fi. The resulting reduced problem 
is an Inverse Poisson problem, of finding a distribution of localized charges 
adhering to these conditions. The charges correspond to cone points. 

The main component needed to apply the theory to mesh generation of planar 
domains is a suitable Inverse Poisson algorithm. An algorithm for creating the 
final discrete mesh is also required. 
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A Relation between geodesic curvature at different metrics 

In this appendix Eq. (4) is derived. Eq. (4) is a known relation in conformal 
geometry, however sign and direction convensions vary, so a derivation is given 
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for completeness. Let (x 1 ,^ 2 ) be a coordinate system, the metric compo- 
nents, and g = det (gij)- For a curve a parametrized by the length parameter 
a (s) = (a 1 (s) , a 2 (s)), n g is given by [16] 

da 1 ( d 2 a % ^ • da j da k \ , . , . 



where en satisfies en = e 2 2 = 0, e i2 = —£21 = y/9- ^jk are the Christoffel 
symbols, related to the metric by the formulas 

r 1 _ l „ki( d 9ik d gij t dg kj \ 

L " ~2 9 [dxi dx k + d&)' 1 J 

where g kl is the k, I entry of the inverse of the matrix (gij). 

Denote by g^ the "standard" metric induced by the embedding in three di- 
mensional space, by g ik the metric of the 0-manifold, (see Definition 3), and 
the corresponding Christoffel symbols by T^^T^ 1 respectively. Substituting 
gij = e 2<i> gij (see Definition 3) into Eq. (A. 2), the equation for T^ 1 reads 

V = r iJ . , + 0, J -(j{-0, fc ^%- + 0, i (jJ , (A.3) 

where comas denote differentiation: 0,j = d<f)/dx l . Eq. (A.3) can now be used 
to derive a relation between K g and K g , the geodesic curvatures in the two 
metrics, at a point p on the curve a. The derivation is simpler if a normal 
coordinate system is chosen, for which g^ (p) = 5ij, and (p) = (this can 
always be done, see [16]). Noting that y/g = e^y/g, ds = e^ds, the equation 



for Kg reads 



da 1 ( d 2 a % ~ „■ da j da k 



kg ~ ~ EU ds [ ds 2 + T > k ds ds ) 

= e ~ % -!7 [l? + (*•' *' ~ + *" *) dTdT j 

,( n . da k ( da 1 da 1 \ da 1 „■„ , / da J da^\\ 
= e * [ Kg + 20, fc — -T—r-eu - —eu5 im <j), m - — 

\ ds \ ds ds J ds \ ds ds J J 

= e~+ (k 9 - d n ct>) . 

The last equality holds because ^f-^f~£ji = , as can be directly verified, be- 
cause ^-^f- = 1 in arc length parametrization. The convension that (T a , N a ) 
form a right-hand system has been used. 
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B Triangular Meshes 



In this appendix the changes in the definitions and results in sections 3,4, and 
5 for the case of triangular meshes are outlined. 

Definition 9 (triangle cross field) Define an equivalence ~ of vectors in 
R 2 : for 2 vectors v 1; v 2 G M, 2 , v x ~ v 2 if and only if v l7 v 2 are form an angle 
of mi/3, with n e Z. A cross is an element ofR 2 / ~. 

The cross-field and 0-manifold are the same as for the quadrilateral case. 
The conditions for the existence of a triangle cross-field are: 
Condition 1 for triangular meshes reads: 
Condition 1. obeys the equation 

V 2 Hr)=K + \ £ k4% 

6 i=l..N 

the Poisson equation with point sources (delta functions) with ki e Z, 

h > -6. 

Condition 2. Is the same as in the quadrilateral case. 

Condition 3. Eq. (36) for the triangular case reads 

f 71 

/ d n <j)ds = n — h C . 



Condition 4. Eq. (37) of becomes 

/ d n (j)ds = 6 a2 - 6 ai +/ K g ds + n-. 

J a± J ai 3 

Fig. B.l shows selected triangle cross-field geodesies around a singularity (with 
(f>L — 0), spaced at manifold-distances of As = 1 from each other, for different 
singularity strengths. 



C Higher-genus surfaces, and a proof of Theorem 4 

This appendix includes a proof of Theorem 4, but first, the case of surfaces 
with genus higher than zero is briefly disscussed. Basic notions of algebraic- 
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(i) (ii) 




Fig. B.l. Triangle cross field geodesies around a singularity. The geodesies are spaced 
at unit manifold-distances, (i) A k = — 1 singularity, (ii) A k = 1 singularity. 

topology are used, see e.g. [22]. 

As the proof of Theorem 4 presented below shows, for a cross- field to exist, 
the angle change in the parallel-transport of a vector along a loop based at ao 
must be a multiple of tt/2. Let D' be the surface D with holes are cut around 
the cone points of D. Let ao be a base point, as in Theorem 4. Let {a{\ be a set 
of loops starting at ao, each encircling a single boundary of D' (including the 
holes cut around the cone points). The angle change due to parallel transport 
around these loops is a multiple of 7r/2, see the proof of the theorem below. 
Complete the set {a{} to a homotopy basis [22] of D f , by adding a set of 
loops {/3i} based at a . The additional condition is that the angle change due 
to parallel transport around a loop of {$} is a multiple of n/2. Any loop 
based at a is homotopic to a composition of loops from the homotopy basis. 
Since the manifold D' is flat with the conformal metric, parallel transport is 
preserved by the homotopy, hence the change in the angle along any loop is 
a multiple of 7r/2 and the cross-field is well-defined. The second part of the 
proof below remains unchanged. 

We now turn to the proof of Theorem 4: 

Proof. We first prove that V (b) for some b G D\P is independent of a. Let 
«i,a 2 be two curves from a to 6. Denote T = (a Q ). We need to show 
that the parallel transport of To to b gives vectors that belong to the same 
cross, i.e. 

l(PT a2 T , Pf ai T ) =k 1X - (C.l) 
for some fc 6 Z. Using Eq. (4), (5), Eq. (C.l) becomes 

k \ = L y 9 - f£) ds - L p - f£) ds = £ { Kg ~ i£) ds (c - 2) 

where a = [a^ajf]. (Note that Eq. (C.2) is equivalent to the requirement 
that the cross in a be parallel-translated to itself along a.) We now prove 
Eq. (C.2). The region enclosed by a can contain singularities and boundary 
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Fig. C.l. Proving that the cross-field is well defined. T r is a boundary curve, p is a 
cone point, ps is a cone point on the boundary. 

curves. Surround them by additional curves, as shown in Fig. C.l. Denote the 
region between the a curve and the u and ip curves as S (here the assumption 
that D has genus zero enters). Since V 2 = K, Eq. (14), everywhere in S, 
and are finite on the boundary of S (due to 0- manifold definition, (hi)), 
Green's theorem, Eq. (8), can be applied in this case, giving 

-l s Kda= -L v2 * ia= iiL is+ £ £^ ds+ £ Ltn is{c - 3) 

Jb Jb J a Oil AT Jcoi (Jfl Ar Jipj Uil 

where N u are the number of points in P, and inner-curves enclosed by 
a. Consider a single loop ip r around an inner curve T r , 1 < r < N^. ip r is 
composed of curves r]i around points of J U P on the boundary, and curves 
Q between junctions, i.e. ip r = rji, (2, f)2, ■■■], see Fig. C.l. The total flux 
through ip r is given by 

L d ^ ds = ? (£ ^ ds + £ §H = ? <S +k 4 + £ K ° ds ) = 

= E(( 7r -^) + ^ + £ Kgds ) 

= E(^+£s^ + ^|) (C.4) 

with kj i: k^ r E Z. The second equality uses Conditions 2,3. The flux through 
uji is given by Condition 1: 

£>-4 
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a = [jui , t?i , ^2] 



Fig. C.2. Proving that the cross-field is aligned with the boundary. 
After substituting Eq. (C.4 ),(C5), Eq. (C.3) becomes 



and using Gauss-Bonnet theorem for a domain with more then one boundary 
component we find 



for some k G Z. This proves Eq. (C.2). 

We now turn to prove the property (ii) of a cross-field, its boundary alignment. 
Let c be a boundary point, c G r r \ (JUP), 1 < r < N^. As was shown 
above, the cross at c, V(c), is well-defined. It is left to show that Tr^ (c) 
G V(c). Denote T r = Tr r (e r ). The assumptions of the theorem, together 
with Condition 4, assure that T r G V (e r ). Let /3 be the section of r r between 
e r and c. Define a curve a from e r to c composed of smooth curves r\ n around 
the points of J U P in the image of /3, and \x m curves between junctions, see 
Fig. C.2. The curves i] n avoid the sigularities that the function <fi might have at 
points in JUP. Let 9j n be the junction angle at the point of JUP "bypassed" 
by the curve rj n (zero for points that are not juction points). Then, if r\ n follows 
T r closely, the total turn of r\ n is equal to 9 Jn : 



as can be formally verified by applying the Gauss-Bonet theorem to the area 
between rji and r r . Note furthermore that 9j n + n = 9 irin , the inner angle at 
that point. 




i=l..N w i=l..N^ 
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[n g - d n 4>) ds 



The variation of Tr> (e r ) along a is 
Z(PT a T r ,T r ,.(c)) = / k g d~s= I 

IT-] J Um ,' ^ "J Tin 



with k e Z. Conditions 2 and 3 were applied for the /x m and % curves, re- 
spectively. Thus T ri (c) G V (C), which completes the proof of the boundary- 
properties of the cross-field 1/. ■ 
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